///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/**
*	Contains code for 4x4 matrices.
*	\file		LABORMatrix4x4.h
*	\author		Pierre Terdiman
*	\date		April, 4, 2000
*/
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Include Guard
#ifndef __LABORMATRIX4X4_H__
#define __LABORMATRIX4X4_H__

#include "LaborMemoryMacros.h"
#include "LaborMatrix3x3.h"
#include "../LaborCoreDef.h"
#include "LaborFPU.h"
#include "LaborMathDecl.h"
#include "LaborPoint4.h"
NAMESPACE_LABOR3D_BEGIN
// Forward declarations
class CLPlane;

#define	MATRIX4X4_EPSILON		(1.0e-7f)

class  LABORCORE_API CLMatrix4x4
{
	//				void	LUBackwardSubstitution( sdword *indx, float* b );
	//				void	LUDecomposition( sdword* indx, float* d );
public:
	static const CLMatrix4x4 ZERO;
	static const CLMatrix4x4 IDENTITY;
	static const CLMatrix4x4 CLIPSPACE2DTOIMAGESPACE;

public:
	//! Empty constructor.
	inline						CLMatrix4x4()									{}
	///< Constructor rotate about origin by point
	CLMatrix4x4( CLPoint3 point );
	//! Constructor from 16 values
	inline						CLMatrix4x4(	float m00, float m01, float m02, float m03,
		float m10, float m11, float m12, float m13,
		float m20, float m21, float m22, float m23,
		float m30, float m31, float m32, float m33)
	{
		m[0][0] = m00;	m[0][1] = m01;	m[0][2] = m02;	m[0][3] = m03;
		m[1][0] = m10;	m[1][1] = m11;	m[1][2] = m12;	m[1][3] = m13;
		m[2][0] = m20;	m[2][1] = m21;	m[2][2] = m22;	m[2][3] = m23;
		m[3][0] = m30;	m[3][1] = m31;	m[3][2] = m32;	m[3][3] = m33;
	}
	//! Copy constructor
	inline						CLMatrix4x4(const CLMatrix4x4& mat)				{ CopyMemory(m, &mat.m, 16*sizeof(float));	}
	inline CLMatrix4x4(float _m[16])
	{
		memcpy((char*)data,(char*)_m,sizeof(float)*16);
	}
	//! Destructor.
	inline						~CLMatrix4x4()								{}

	//! Assign values (rotation only)
	inline	CLMatrix4x4&			set(	float m00, float m01, float m02,
		float m10, float m11, float m12,
		float m20, float m21, float m22)
	{
		m[0][0] = m00;	m[0][1] = m01;	m[0][2] = m02;
		m[1][0] = m10;	m[1][1] = m11;	m[1][2] = m12;
		m[2][0] = m20;	m[2][1] = m21;	m[2][2] = m22;
		return *this;
	}
	//! Assign values
	inline	CLMatrix4x4&			set(	float m00, float m01, float m02, float m03,
		float m10, float m11, float m12, float m13,
		float m20, float m21, float m22, float m23,
		float m30, float m31, float m32, float m33)
	{
		m[0][0] = m00;	m[0][1] = m01;	m[0][2] = m02;	m[0][3] = m03;
		m[1][0] = m10;	m[1][1] = m11;	m[1][2] = m12;	m[1][3] = m13;
		m[2][0] = m20;	m[2][1] = m21;	m[2][2] = m22;	m[2][3] = m23;
		m[3][0] = m30;	m[3][1] = m31;	m[3][2] = m32;	m[3][3] = m33;
		return *this;
	}

	//! get indenty matrix
	static	const CLMatrix4x4&	getIndentity();

	//! Copy from a CLMatrix4x4
	inline	void				copy(const CLMatrix4x4& source)				{ CopyMemory(m, source.m, 16*sizeof(float));	}

	// Row-column access
	//! Returns a row.
	inline	void				getRow(const unsigned int r, CLPoint4& p)	const	{ p.x=m[r][0];	p.y=m[r][1];	p.z=m[r][2];	p.w=m[r][3];		}
	//! Returns a row.
	inline	void				getRow(const unsigned int r, CLPoint3& p)	const	{ p.x=m[r][0];	p.y=m[r][1];	p.z=m[r][2];						}
	//! Returns a row.
	inline	const CLPoint4&		getRow(const unsigned int r)				const	{ return *(const CLPoint4*)&m[r][0];									}
	//! Returns a row.
	inline	CLPoint4&			getRow(const unsigned int r)						{ return *(CLPoint4*)&m[r][0];										}
	//! Sets a row.
	inline	void				setRow(const unsigned int r, const CLPoint4& p)		{ m[r][0]=p.x;	m[r][1]=p.y;	m[r][2]=p.z;	m[r][3]=p.w;		}
	//! Sets a row.
	inline	void				setRow(const unsigned int r, const CLPoint3& p)		{ m[r][0]=p.x;	m[r][1]=p.y;	m[r][2]=p.z;	m[r][3]= (r!=3) ? 0.0f : 1.0f;		}
	//! Returns a column.
	inline	void				getCol(const unsigned int c, CLPoint4& p)	const	{ p.x=m[0][c];	p.y=m[1][c];	p.z=m[2][c];	p.w=m[3][c];		}
	//! Returns a column.
	inline	void				getCol(const unsigned int c, CLPoint3& p)	const	{ p.x=m[0][c];	p.y=m[1][c];	p.z=m[2][c];						}
	//! Sets a column.
	inline	void				setCol(const unsigned int c, const CLPoint4& p)		{ m[0][c]=p.x;	m[1][c]=p.y;	m[2][c]=p.z;	m[3][c]=p.w;		}
	//! Sets a column.
	inline	void				setCol(const unsigned int c, const CLPoint3& p)		{ m[0][c]=p.x;	m[1][c]=p.y;	m[2][c]=p.z;	m[3][c]= (c!=3) ? 0.0f : 1.0f;	}

	// Translation
	//! Returns the translation part of the matrix.
	inline	const CLPoint4&		getTrans()							const	{ return getRow(3);								}
	//! Gets the translation part of the matrix
	inline	void				getTrans(CLPoint3& p)					const	{ p.x=m[3][0];	p.y=m[3][1];	p.z=m[3][2];	}
	//! Sets the translation part of the matrix, from a Point.
	inline	void				setTrans(const CLPoint3& p)					{ m[3][0]=p.x;	m[3][1]=p.y;	m[3][2]=p.z;	}
	//! Sets the translation part of the matrix, from a CLPoint4.
	inline	void				setTrans(const CLPoint4& p)					{ m[3][0]=p.x;	m[3][1]=p.y;	m[3][2]=p.z;	m[3][3]=p.w;	}
	//! Sets the translation part of the matrix, from floats.
	inline	void				setTrans(float tx, float ty, float tz)		{ m[3][0]=tx;	m[3][1]=ty;		m[3][2]=tz;		}

	// Scale
	//! Sets the scale from a Point. The point is put on the diagonal.
	inline	void				setScale(const CLPoint3& p)					{ m[0][0]=p.x;	m[1][1]=p.y;	m[2][2]=p.z;	}
	//! Sets the scale from floats. Values are put on the diagonal.
	inline	void				setScale(float sx, float sy, float sz)		{ m[0][0]=sx;	m[1][1]=sy;		m[2][2]=sz;		}
	//! Scales from a Point. Each row is multiplied by a component.
	void				scale(const CLPoint3& p)
	{
		m[0][0] *= p.x;	m[1][0] *= p.y;	m[2][0] *= p.z;
		m[0][1] *= p.x;	m[1][1] *= p.y;	m[2][1] *= p.z;
		m[0][2] *= p.x;	m[1][2] *= p.y;	m[2][2] *= p.z;
	}
	//! Scales from floats. Each row is multiplied by a value.
	void				scale(float sx, float sy, float sz)
	{
		m[0][0] *= sx;	m[1][0] *= sy;	m[2][0] *= sz;
		m[0][1] *= sx;	m[1][1] *= sy;	m[2][1] *= sz;
		m[0][2] *= sx;	m[1][2] *= sy;	m[2][2] *= sz;
	}
	void                noScale();

	//! Computes the trace. The trace is the sum of the 4 diagonal components.
	inline	float				trace()			const			{ return m[0][0] + m[1][1] + m[2][2] + m[3][3];			}
	//! Computes the trace of the upper 3x3 matrix.
	inline	float				trace3x3()		const			{ return m[0][0] + m[1][1] + m[2][2];					}
	//! Clears the matrix.
	inline	void				zero()							{ ZeroMemory(&m,  sizeof(m));							}
	//! Sets the identity matrix.
	inline	void				identity()						{ zero(); m[0][0] = m[1][1] = m[2][2] = m[3][3] = 1.0f;	}
	//! Checks for identity
	inline	bool				isIdentity()	const
	{
		if(IR(m[0][0])!=IEEE_1_0)	return false;
		if(IR(m[0][1])!=0)			return false;
		if(IR(m[0][2])!=0)			return false;
		if(IR(m[0][3])!=0)			return false;

		if(IR(m[1][0])!=0)			return false;
		if(IR(m[1][1])!=IEEE_1_0)	return false;
		if(IR(m[1][2])!=0)			return false;
		if(IR(m[1][3])!=0)			return false;

		if(IR(m[2][0])!=0)			return false;
		if(IR(m[2][1])!=0)			return false;
		if(IR(m[2][2])!=IEEE_1_0)	return false;
		if(IR(m[2][3])!=0)			return false;

		if(IR(m[3][0])!=0)			return false;
		if(IR(m[3][1])!=0)			return false;
		if(IR(m[3][2])!=0)			return false;
		if(IR(m[3][3])!=IEEE_1_0)	return false;
		return true;
	}

	//! Checks matrix validity
	inline	bool				isValid()		const
	{
		for(unsigned int j=0;j<4;j++)
		{
			for(unsigned int i=0;i<4;i++)
			{
				if(!isValidFloat(m[j][i]))	return false;
			}
		}
		return true;
	}

	//! Sets a rotation matrix around the X axis.
	void				rotX(float angle)	{ float Cos = cosf(angle), Sin = sinf(angle); identity(); m[1][1] = m[2][2] = Cos; m[2][1] = -Sin;	m[1][2] = Sin;	}
	//! Sets a rotation matrix around the Y axis.
	void				rotY(float angle)	{ float Cos = cosf(angle), Sin = sinf(angle); identity(); m[0][0] = m[2][2] = Cos; m[2][0] = Sin;	m[0][2] = -Sin;	}
	//! Sets a rotation matrix around the Z axis.
	void				rotZ(float angle)	{ float Cos = cosf(angle), Sin = sinf(angle); identity(); m[0][0] = m[1][1] = Cos; m[1][0] = -Sin;	m[0][1] = Sin;	}

	//! Makes a rotation matrix about an arbitrary axis
	CLMatrix4x4&			rot(float angle, const CLPoint3& axis );

	CLMatrix4x4&			buildReflectionMatrix( CLPlane& p );
	CLMatrix4x4&			buildReflectionMatrix( const CLPlane& p );


	void				blend( float fFactor, const CLMatrix4x4 &matrix ){
		m[0][0] += fFactor*matrix.m[0][0]; m[0][1] += fFactor*matrix.m[0][1]; m[0][2] += fFactor*matrix.m[0][2]; m[0][3] += fFactor*matrix.m[0][3];
		m[1][0] += fFactor*matrix.m[1][0]; m[1][1] += fFactor*matrix.m[1][1]; m[0][2] += fFactor*matrix.m[1][2]; m[1][3] += fFactor*matrix.m[0][3];
		m[2][0] += fFactor*matrix.m[2][0]; m[2][1] += fFactor*matrix.m[2][1]; m[0][2] += fFactor*matrix.m[2][2]; m[2][3] += fFactor*matrix.m[0][3];
		m[3][0] += fFactor*matrix.m[3][0]; m[3][1] += fFactor*matrix.m[3][1]; m[0][2] += fFactor*matrix.m[3][2]; m[3][3] += fFactor*matrix.m[0][3];
	}

	/// build matix left hand
	const CLMatrix4x4	&lookAtLH( const CLPoint3* pvPos, const CLPoint3* pvAt , const CLPoint3* pvUp ){
		CLPoint3 vDir,vRight,vUp;
		vDir = *pvAt - *pvPos;
		vDir.normalize();

		CLPoint3 vPos = *pvPos;

		vRight = pvUp->cross( vDir );
		vRight.normalize();
		vUp = vDir.cross( vRight );

		//CLMatrix4x4 matrixResult;

		m[0][0] = vRight.x;    m[0][1] = vUp.x;    m[0][2] = vDir.x;
		m[1][0] = vRight.y;    m[1][1] = vUp.y;    m[1][2] = vDir.y;
		m[2][0] = vRight.z;    m[2][1] = vUp.z;    m[2][2] = vDir.z;

		m[3][0]	= -vPos.dot( vRight );
		m[3][1]	= -vPos.dot( vUp );
		m[3][2]	= -vPos.dot( vDir );

		m[0][3]	= 0.f;
		m[1][3]	= 0.f;
		m[2][3]	= 0.f;
		m[3][3]	= 1.f;

		return *this;
	}

	//! Transposes the matrix.
	void				transpose()
	{
		IR(m[1][0]) ^= IR(m[0][1]);		IR(m[0][1]) ^= IR(m[1][0]);		IR(m[1][0]) ^= IR(m[0][1]);
		IR(m[2][0]) ^= IR(m[0][2]);		IR(m[0][2]) ^= IR(m[2][0]);		IR(m[2][0]) ^= IR(m[0][2]);
		IR(m[3][0]) ^= IR(m[0][3]);		IR(m[0][3]) ^= IR(m[3][0]);		IR(m[3][0]) ^= IR(m[0][3]);
		IR(m[1][2]) ^= IR(m[2][1]);		IR(m[2][1]) ^= IR(m[1][2]);		IR(m[1][2]) ^= IR(m[2][1]);
		IR(m[1][3]) ^= IR(m[3][1]);		IR(m[3][1]) ^= IR(m[1][3]);		IR(m[1][3]) ^= IR(m[3][1]);
		IR(m[2][3]) ^= IR(m[3][2]);		IR(m[3][2]) ^= IR(m[2][3]);		IR(m[2][3]) ^= IR(m[3][2]);
	}

	//! Computes a cofactor. Used for matrix inversion.
	float				coFactor(unsigned int row, unsigned int col)	const;
	//! Computes the determinant of the matrix.
	float				determinant()	const;
	//! Inverts the matrix. Determinant must be different from zero, else matrix can't be inverted.
	CLMatrix4x4&			invert();
	CLMatrix4x4&			affineInvert();
	//				Matrix&	ComputeAxisMatrix(CLPoint3& axis, float angle);

	// Cast operators
	//! Casts a CLMatrix4x4 to a Matrix3x3.
	inline	operator			CLMatrix3x3()	const
	{
		return CLMatrix3x3(
			m[0][0],	m[0][1],	m[0][2],
			m[1][0],	m[1][1],	m[1][2],
			m[2][0],	m[2][1],	m[2][2]);
	}
	//! Casts a CLMatrix4x4 to a CLQuat.
	operator			CLQuat()	const;
	// Arithmetic operators
	//! Operator for CLMatrix4x4 Plus = CLMatrix4x4 + CLMatrix4x4;
	inline	CLMatrix4x4			operator+(const CLMatrix4x4& mat)	const
	{
		return CLMatrix4x4(
			m[0][0]+mat.m[0][0], m[0][1]+mat.m[0][1], m[0][2]+mat.m[0][2], m[0][3]+mat.m[0][3], 
			m[1][0]+mat.m[1][0], m[1][1]+mat.m[1][1], m[1][2]+mat.m[1][2], m[1][3]+mat.m[1][3], 
			m[2][0]+mat.m[2][0], m[2][1]+mat.m[2][1], m[2][2]+mat.m[2][2], m[2][3]+mat.m[2][3], 
			m[3][0]+mat.m[3][0], m[3][1]+mat.m[3][1], m[3][2]+mat.m[3][2], m[3][3]+mat.m[3][3]);
	}

	//! Operator for CLMatrix4x4 Minus = CLMatrix4x4 - CLMatrix4x4;
	inline	CLMatrix4x4			operator-(const CLMatrix4x4& mat)	const
	{
		return CLMatrix4x4(
			m[0][0]-mat.m[0][0], m[0][1]-mat.m[0][1], m[0][2]-mat.m[0][2], m[0][3]-mat.m[0][3], 
			m[1][0]-mat.m[1][0], m[1][1]-mat.m[1][1], m[1][2]-mat.m[1][2], m[1][3]-mat.m[1][3], 
			m[2][0]-mat.m[2][0], m[2][1]-mat.m[2][1], m[2][2]-mat.m[2][2], m[2][3]-mat.m[2][3], 
			m[3][0]-mat.m[3][0], m[3][1]-mat.m[3][1], m[3][2]-mat.m[3][2], m[3][3]-mat.m[3][3]);
	}

	//! Operator for CLMatrix4x4 Mul = CLMatrix4x4 * CLMatrix4x4;
	inline	CLMatrix4x4			operator*(const CLMatrix4x4& mat)	const
	{
		return CLMatrix4x4(
			m[0][0]*mat.m[0][0] + m[0][1]*mat.m[1][0] + m[0][2]*mat.m[2][0] + m[0][3]*mat.m[3][0],
			m[0][0]*mat.m[0][1] + m[0][1]*mat.m[1][1] + m[0][2]*mat.m[2][1] + m[0][3]*mat.m[3][1],
			m[0][0]*mat.m[0][2] + m[0][1]*mat.m[1][2] + m[0][2]*mat.m[2][2] + m[0][3]*mat.m[3][2],
			m[0][0]*mat.m[0][3] + m[0][1]*mat.m[1][3] + m[0][2]*mat.m[2][3] + m[0][3]*mat.m[3][3],

			m[1][0]*mat.m[0][0] + m[1][1]*mat.m[1][0] + m[1][2]*mat.m[2][0] + m[1][3]*mat.m[3][0],
			m[1][0]*mat.m[0][1] + m[1][1]*mat.m[1][1] + m[1][2]*mat.m[2][1] + m[1][3]*mat.m[3][1],
			m[1][0]*mat.m[0][2] + m[1][1]*mat.m[1][2] + m[1][2]*mat.m[2][2] + m[1][3]*mat.m[3][2],
			m[1][0]*mat.m[0][3] + m[1][1]*mat.m[1][3] + m[1][2]*mat.m[2][3] + m[1][3]*mat.m[3][3],

			m[2][0]*mat.m[0][0] + m[2][1]*mat.m[1][0] + m[2][2]*mat.m[2][0] + m[2][3]*mat.m[3][0],
			m[2][0]*mat.m[0][1] + m[2][1]*mat.m[1][1] + m[2][2]*mat.m[2][1] + m[2][3]*mat.m[3][1],
			m[2][0]*mat.m[0][2] + m[2][1]*mat.m[1][2] + m[2][2]*mat.m[2][2] + m[2][3]*mat.m[3][2],
			m[2][0]*mat.m[0][3] + m[2][1]*mat.m[1][3] + m[2][2]*mat.m[2][3] + m[2][3]*mat.m[3][3],

			m[3][0]*mat.m[0][0] + m[3][1]*mat.m[1][0] + m[3][2]*mat.m[2][0] + m[3][3]*mat.m[3][0],
			m[3][0]*mat.m[0][1] + m[3][1]*mat.m[1][1] + m[3][2]*mat.m[2][1] + m[3][3]*mat.m[3][1],
			m[3][0]*mat.m[0][2] + m[3][1]*mat.m[1][2] + m[3][2]*mat.m[2][2] + m[3][3]*mat.m[3][2],
			m[3][0]*mat.m[0][3] + m[3][1]*mat.m[1][3] + m[3][2]*mat.m[2][3] + m[3][3]*mat.m[3][3]);
	}

	//! Operator for CLPoint4 Mul = CLMatrix4x4 * CLPoint4;
	inline	CLPoint4				operator*(const CLPoint4& v)		const	{ return CLPoint4(getRow(0)|v, getRow(1)|v, getRow(2)|v, getRow(3)|v); }

	//! Operator for Point Mul = CLMatrix4x4 * Point;
	inline	CLPoint3				operator*(const CLPoint3& v)		const
	{
		return CLPoint3(	m[0][0]*v.x + m[0][1]*v.y + m[0][2]*v.z + m[0][3],
			m[1][0]*v.x + m[1][1]*v.y + m[1][2]*v.z + m[1][3],
			m[2][0]*v.x + m[2][1]*v.y + m[2][2]*v.z + m[2][3]	);
	}

	//! Operator for Matrix4x4 Scale = CLMatrix4x4 * float;
	inline	CLMatrix4x4			operator*(float s)				const
	{
		return CLMatrix4x4(
			m[0][0]*s,	m[0][1]*s,	m[0][2]*s,	m[0][3]*s,
			m[1][0]*s,	m[1][1]*s,	m[1][2]*s,	m[1][3]*s,
			m[2][0]*s,	m[2][1]*s,	m[2][2]*s,	m[2][3]*s,
			m[3][0]*s,	m[3][1]*s,	m[3][2]*s,	m[3][3]*s);
	}

	//! Operator for CLMatrix4x4 Scale = float * Matrix4x4;
	inline	friend CLMatrix4x4	operator*(float s, const CLMatrix4x4& mat)
	{
		return CLMatrix4x4(
			s*mat.m[0][0],	s*mat.m[0][1],	s*mat.m[0][2],	s*mat.m[0][3],
			s*mat.m[1][0],	s*mat.m[1][1],	s*mat.m[1][2],	s*mat.m[1][3],
			s*mat.m[2][0],	s*mat.m[2][1],	s*mat.m[2][2],	s*mat.m[2][3],
			s*mat.m[3][0],	s*mat.m[3][1],	s*mat.m[3][2],	s*mat.m[3][3]);
	}

	//! Operator for Matrix4x4 Div = Matrix4x4 / float;
	inline	CLMatrix4x4			operator/(float s)				const
	{
		if(s) s = 1.0f / s;

		return CLMatrix4x4(
			m[0][0]*s,	m[0][1]*s,	m[0][2]*s,	m[0][3]*s,
			m[1][0]*s,	m[1][1]*s,	m[1][2]*s,	m[1][3]*s,
			m[2][0]*s,	m[2][1]*s,	m[2][2]*s,	m[2][3]*s,
			m[3][0]*s,	m[3][1]*s,	m[3][2]*s,	m[3][3]*s);
	}

	//! Operator for Matrix4x4 Div = float / Matrix4x4;
	inline	friend CLMatrix4x4	operator/(float s, const CLMatrix4x4& mat)
	{
		return CLMatrix4x4(
			s/mat.m[0][0],	s/mat.m[0][1],	s/mat.m[0][2],	s/mat.m[0][3],
			s/mat.m[1][0],	s/mat.m[1][1],	s/mat.m[1][2],	s/mat.m[1][3],
			s/mat.m[2][0],	s/mat.m[2][1],	s/mat.m[2][2],	s/mat.m[2][3],
			s/mat.m[3][0],	s/mat.m[3][1],	s/mat.m[3][2],	s/mat.m[3][3]);
	}

	//! Operator for Matrix4x4 += Matrix4x4;
	inline	CLMatrix4x4&			operator+=(const CLMatrix4x4& mat)
	{
		m[0][0]+=mat.m[0][0];	m[0][1]+=mat.m[0][1];	m[0][2]+=mat.m[0][2];	m[0][3]+=mat.m[0][3];
		m[1][0]+=mat.m[1][0];	m[1][1]+=mat.m[1][1];	m[1][2]+=mat.m[1][2];	m[1][3]+=mat.m[1][3];
		m[2][0]+=mat.m[2][0];	m[2][1]+=mat.m[2][1];	m[2][2]+=mat.m[2][2];	m[2][3]+=mat.m[2][3];
		m[3][0]+=mat.m[3][0];	m[3][1]+=mat.m[3][1];	m[3][2]+=mat.m[3][2];	m[3][3]+=mat.m[3][3];
		return	*this;
	}

	//! Operator for Matrix4x4 -= Matrix4x4;
	inline	CLMatrix4x4&			operator-=(const CLMatrix4x4& mat)
	{
		m[0][0]-=mat.m[0][0];	m[0][1]-=mat.m[0][1];	m[0][2]-=mat.m[0][2];	m[0][3]-=mat.m[0][3];
		m[1][0]-=mat.m[1][0];	m[1][1]-=mat.m[1][1];	m[1][2]-=mat.m[1][2];	m[1][3]-=mat.m[1][3];
		m[2][0]-=mat.m[2][0];	m[2][1]-=mat.m[2][1];	m[2][2]-=mat.m[2][2];	m[2][3]-=mat.m[2][3];
		m[3][0]-=mat.m[3][0];	m[3][1]-=mat.m[3][1];	m[3][2]-=mat.m[3][2];	m[3][3]-=mat.m[3][3];
		return	*this;
	}

	//! Operator for Matrix4x4 *= Matrix4x4;
	CLMatrix4x4&			operator*=(const CLMatrix4x4& mat)
	{
		CLPoint4 TempRow;

		getRow(0, TempRow);
		m[0][0] = TempRow.x*mat.m[0][0] + TempRow.y*mat.m[1][0] + TempRow.z*mat.m[2][0] + TempRow.w*mat.m[3][0];
		m[0][1] = TempRow.x*mat.m[0][1] + TempRow.y*mat.m[1][1] + TempRow.z*mat.m[2][1] + TempRow.w*mat.m[3][1];
		m[0][2] = TempRow.x*mat.m[0][2] + TempRow.y*mat.m[1][2] + TempRow.z*mat.m[2][2] + TempRow.w*mat.m[3][2];
		m[0][3] = TempRow.x*mat.m[0][3] + TempRow.y*mat.m[1][3] + TempRow.z*mat.m[2][3] + TempRow.w*mat.m[3][3];

		getRow(1, TempRow);
		m[1][0] = TempRow.x*mat.m[0][0] + TempRow.y*mat.m[1][0] + TempRow.z*mat.m[2][0] + TempRow.w*mat.m[3][0];
		m[1][1] = TempRow.x*mat.m[0][1] + TempRow.y*mat.m[1][1] + TempRow.z*mat.m[2][1] + TempRow.w*mat.m[3][1];
		m[1][2] = TempRow.x*mat.m[0][2] + TempRow.y*mat.m[1][2] + TempRow.z*mat.m[2][2] + TempRow.w*mat.m[3][2];
		m[1][3] = TempRow.x*mat.m[0][3] + TempRow.y*mat.m[1][3] + TempRow.z*mat.m[2][3] + TempRow.w*mat.m[3][3];

		getRow(2, TempRow);
		m[2][0] = TempRow.x*mat.m[0][0] + TempRow.y*mat.m[1][0] + TempRow.z*mat.m[2][0] + TempRow.w*mat.m[3][0];
		m[2][1] = TempRow.x*mat.m[0][1] + TempRow.y*mat.m[1][1] + TempRow.z*mat.m[2][1] + TempRow.w*mat.m[3][1];
		m[2][2] = TempRow.x*mat.m[0][2] + TempRow.y*mat.m[1][2] + TempRow.z*mat.m[2][2] + TempRow.w*mat.m[3][2];
		m[2][3] = TempRow.x*mat.m[0][3] + TempRow.y*mat.m[1][3] + TempRow.z*mat.m[2][3] + TempRow.w*mat.m[3][3];

		getRow(3, TempRow);
		m[3][0] = TempRow.x*mat.m[0][0] + TempRow.y*mat.m[1][0] + TempRow.z*mat.m[2][0] + TempRow.w*mat.m[3][0];
		m[3][1] = TempRow.x*mat.m[0][1] + TempRow.y*mat.m[1][1] + TempRow.z*mat.m[2][1] + TempRow.w*mat.m[3][1];
		m[3][2] = TempRow.x*mat.m[0][2] + TempRow.y*mat.m[1][2] + TempRow.z*mat.m[2][2] + TempRow.w*mat.m[3][2];
		m[3][3] = TempRow.x*mat.m[0][3] + TempRow.y*mat.m[1][3] + TempRow.z*mat.m[2][3] + TempRow.w*mat.m[3][3];

		return	*this;
	}

	//! Operator for Matrix4x4 *= float;
	inline	CLMatrix4x4&		operator*=(float s)
	{
		m[0][0]*=s;	m[0][1]*=s;	m[0][2]*=s;	m[0][3]*=s;
		m[1][0]*=s;	m[1][1]*=s;	m[1][2]*=s;	m[1][3]*=s;
		m[2][0]*=s;	m[2][1]*=s;	m[2][2]*=s;	m[2][3]*=s;
		m[3][0]*=s;	m[3][1]*=s;	m[3][2]*=s;	m[3][3]*=s;
		return	*this;
	}

	//! Operator for Matrix4x4 /= float;
	inline	CLMatrix4x4&		operator/=(float s)
	{
		if(s)  s = 1.0f / s;
		m[0][0]*=s;	m[0][1]*=s;	m[0][2]*=s;	m[0][3]*=s;
		m[1][0]*=s;	m[1][1]*=s;	m[1][2]*=s;	m[1][3]*=s;
		m[2][0]*=s;	m[2][1]*=s;	m[2][2]*=s;	m[2][3]*=s;
		m[3][0]*=s;	m[3][1]*=s;	m[3][2]*=s;	m[3][3]*=s;
		return	*this;
	}
	inline bool            operator==(CLMatrix4x4& mat)
	{
		if(m[0][0] == mat.m[0][0]&&m[0][0] == mat.m[0][1]&&m[0][1] == mat.m[0][1]&&m[0][2] == mat.m[0][2]&&
			m[1][0] == mat.m[1][0]&&m[1][1] == mat.m[1][1]&&m[1][2] == mat.m[1][2]&&m[1][3] == mat.m[1][3]&&
			m[2][0] == mat.m[2][0]&&m[2][1] == mat.m[2][1]&&m[2][2] == mat.m[2][2]&&m[2][3] == mat.m[2][3]&&
			m[3][0] == mat.m[3][0]&&m[3][1] == mat.m[3][1]&&m[3][2] == mat.m[3][2]&&m[3][3] == mat.m[3][3])
			return true;
		else
			return false;
	}
	inline bool            operator==(const CLMatrix4x4& mat)
	{
		if(m[0][0] == mat.m[0][0]&&m[0][0] == mat.m[0][1]&&m[0][1] == mat.m[0][1]&&m[0][2] == mat.m[0][2]&&
			m[1][0] == mat.m[1][0]&&m[1][1] == mat.m[1][1]&&m[1][2] == mat.m[1][2]&&m[1][3] == mat.m[1][3]&&
			m[2][0] == mat.m[2][0]&&m[2][1] == mat.m[2][1]&&m[2][2] == mat.m[2][2]&&m[2][3] == mat.m[2][3]&&
			m[3][0] == mat.m[3][0]&&m[3][1] == mat.m[3][1]&&m[3][2] == mat.m[3][2]&&m[3][3] == mat.m[3][3])
			return true;
		else
			return false;
	}
	inline	const CLPoint4&	operator[](int row)		const	{ return *(const CLPoint4*)&m[row][0];	}
	inline	CLPoint4&			operator[](int row)			{ return *(CLPoint4*)&m[row][0];			}

public:

	union
	{
		float m[4][4];
		float data[16];
#ifdef __cplusplus
		struct
		{
			float m00, m01, m02, m03;
			float m10, m11, m12, m13;
			float m20, m21, m22, m23;
			float m30, m31, m32, m33;
		};
#endif
	};

};

//! Quickly rotates & translates a vector, using the 4x3 part of a 4x4 matrix
inline void transformPoint4x3(CLPoint3& dest, const CLPoint3& source, const CLMatrix4x4& rot)
{
	dest.x = rot.m[3][0] + source.x * rot.m[0][0] + source.y * rot.m[1][0] + source.z * rot.m[2][0];
	dest.y = rot.m[3][1] + source.x * rot.m[0][1] + source.y * rot.m[1][1] + source.z * rot.m[2][1];
	dest.z = rot.m[3][2] + source.x * rot.m[0][2] + source.y * rot.m[1][2] + source.z * rot.m[2][2];
}

//! Quickly rotates a vector, using the 3x3 part of a 4x4 matrix
inline void transformPoint3x3(CLPoint3& dest, const CLPoint3& source, const CLMatrix4x4& rot)
{
	dest.x = source.x * rot.m[0][0] + source.y * rot.m[1][0] + source.z * rot.m[2][0];
	dest.y = source.x * rot.m[0][1] + source.y * rot.m[1][1] + source.z * rot.m[2][1];
	dest.z = source.x * rot.m[0][2] + source.y * rot.m[1][2] + source.z * rot.m[2][2];
}

void invertPRMatrix(CLMatrix4x4& dest, const CLMatrix4x4& src);



//======================================================================///

//        2*zn/w  0       0              0
//	      0       2*zn/h  0              0
//	      0       0       zf/(zn-zf)    -1
//	      0       0       zn*zf/(zn-zf)  0

//======================================================================///
inline void Matrix4x4PerspectiveRH(CLMatrix4x4& mOut, float w, float h, float zn, float zf)
{
	mOut.m01 = mOut.m02 = mOut.m03 = mOut.m10 = mOut.m12 = mOut.m13 = mOut.m20 = mOut.m21 = mOut.m30 = mOut.m31 = mOut.m33 = 0.0f;
	mOut.m00 = 2.0f * zn / w;
	mOut.m11 = 2.0f * zn / h;
	mOut.m22 = zf / (zn - zf);
	mOut.m32 = zn * zf / (zn - zf);
	mOut.m23 = -1.0f;
}

//======================================================================///

//        2*zn/w  0       0              0
//	      0       2*zn/h  0              0
//	      0       0       zf/(zf-zn)     1
//	      0       0       zn*zf/(zn-zf)  0

//======================================================================///

inline void Matrix4x4PerspectiveLH(CLMatrix4x4& mOut, float w, float h, float zn, float zf)
{
	mOut.m01 = mOut.m02 = mOut.m03 = mOut.m10 = mOut.m12 = mOut.m13 = mOut.m20 = mOut.m21 = mOut.m30 = mOut.m31 = mOut.m33 = 0.0f;
	mOut.m00 = 2.0f * zn / w;
	mOut.m11 = 2.0f * zn / h;
	mOut.m22 = zf / (zf - zn);
	mOut.m32 = zn * zf / (zn - zf);
	mOut.m23 = 1.0f;
}

//======================================================================///

//        xScale     0          0              0
//	      0        yScale       0              0
//	      0        0        zf/(zn-zf)        -1
//	      0        0        zn*zf/(zn-zf)      0
//      where:
//      yScale = cot(fovY/2)
        
//      	xScale = yScale / aspect ratio
//======================================================================///
inline void Matrix4x4PerspectiveFovRH(CLMatrix4x4& mOut, float fovy, float aspect, float zn, float zf)
{
	mOut.m01 = mOut.m02 = mOut.m03 = mOut.m10 = mOut.m12 = mOut.m13 = mOut.m20 = mOut.m21 = mOut.m30 = mOut.m31 = mOut.m33 = 0.0f;
	mOut.m11 = 1.0f / (float)tan(fovy * 0.5f);
	mOut.m00 = mOut.m11 / aspect;
	mOut.m22 = zf / (zn - zf);
	mOut.m32 = zn * zf / (zn - zf);
	mOut.m23 = -1.0f;
}

//======================================================================///
//          xScale     0          0               0
//        	0        yScale       0               0
//        	0          0       zf/(zf-zn)         1
//        	0          0       -zn*zf/(zf-zn)     0

//        where:
//        yScale = cot(fovY/2)
//        
//        	xScale = yScale / aspect ratio
//======================================================================///
inline void Matrix4x4PerspectiveFovLH(CLMatrix4x4& mOut, float fovy, float aspect, float zn, float zf)
{
	mOut.m01 = mOut.m02 = mOut.m03 = mOut.m10 = mOut.m12 = mOut.m13 = mOut.m20 = mOut.m21 = mOut.m30 = mOut.m31 = mOut.m33 = 0.0f;
	mOut.m11 = 1.0f / (float)tan(fovy * 0.5f);
	mOut.m00 = mOut.m11 / aspect;
	mOut.m22 = zf / (zf - zn);
	mOut.m32 = zn * zf / (zn - zf);
	mOut.m23 = 1.0f;
}

//======================================================================///

//	        2*zn/(r-l)   0            0                0
//	      	0            2*zn/(t-b)   0                0
//	      	(l+r)/(r-l)  (t+b)/(t-b)  zf/(zn-zf)      -1
//	      	0            0            zn*zf/(zn-zf)    0

//======================================================================///

inline void Matrix4x4PerspectiveOffCenterRH(CLMatrix4x4& mOut, float l, float r, float b, float t, float zn, float zf)
{
	mOut.m01 = mOut.m02 = mOut.m03 = mOut.m10 = mOut.m12 = mOut.m13 = mOut.m30 = mOut.m31 = mOut.m33 = 0.0f;
	mOut.m00 = 2.0f * zn / (r - l);
	mOut.m11 = 2.0f * zn / (t - b);
	mOut.m20 = (l + r) / (l - r);
	mOut.m21 = (t + b) / (b - t);
	mOut.m22 = zf / (zn - zf);
	mOut.m32 = zn * zf / (zn - zf);
	mOut.m23 = -1.0f;
}

//======================================================================///

//	      	    2*zn/(r-l)   0            0              0
//	      		0            2*zn/(t-b)   0              0
//	      		(l+r)/(l-r)  (t+b)/(b-t)  zf/(zf-zn)     1
//	      		0            0            zn*zf/(zn-zf)  0


//======================================================================///

inline void Matrix4x4PerspectiveOffCenterLH(CLMatrix4x4& mOut, float l, float r, float b, float t, float zn, float zf)
{
	mOut.m01 = mOut.m02= mOut.m03 = mOut.m10 = mOut.m12 = mOut.m13 = mOut.m30 = mOut.m31 = mOut.m33 = 0.0f;
	mOut.m00 = 2.0f * zn / (r - l);
	mOut.m11= 2.0f * zn / (t - b);
	mOut.m20 = (l + r) / (l - r);
	mOut.m21 = (t + b) / (b - t);
	mOut.m22 = zf / (zf - zn);
	mOut.m32 = zn * zf / (zn - zf);
	mOut.m23 = 1.0f;
}

//======================================================================///

//	      	2/w  0    0           0
//	      	0    2/h  0           0
//	      	0    0    1/(zn-zf)   0
//	      	0    0    zn/(zn-zf)  1

//======================================================================///

inline void Matrix4x4OrthoRH(CLMatrix4x4& mOut, float w, float h, float zn, float zf)
{
	mOut.m01 = mOut.m02 = mOut.m03 = mOut.m10 = mOut.m12 = mOut.m13 = mOut.m20 = mOut.m21 = mOut.m23 = mOut.m30 = mOut.m31 = 0.0f;
	mOut.m00 = 2.0f / w;
	mOut.m11 = 2.0f / h;
	mOut.m22 = 1.0f / (zn - zf);
	mOut.m32 = zn / (zn - zf);
	mOut.m33 = 1.0f;
}

//======================================================================///

//	      	2/w  0    0           0
//	      	0    2/h  0           0
//	      	0    0    1/(zf-zn)   0
//	      	0    0   -zn/(zf-zn)  1

//======================================================================///

inline void Matrix4x4OrthoLH(CLMatrix4x4& mOut, float w, float h, float zn, float zf)
{
	mOut.m01 = mOut.m02 = mOut.m03 = mOut.m10 = mOut.m12 = mOut.m13 = mOut.m20 = mOut.m21 = mOut.m23 = mOut.m30 = mOut.m31 = 0.0f;
	mOut.m00 = 2.0f / w;
	mOut.m11 = 2.0f / h;
	mOut.m22 = 1.0f / (zf - zn);
	mOut.m32 = zn / (zn - zf);
	mOut.m33 = 1.0f;
}

//======================================================================///

//	      	2/(r-l)      0            0           0
//	      	0            2/(t-b)      0           0
//	      	0            0            1/(zn-zf)   0
//	      	(l+r)/(l-r)  (t+b)/(b-t)  zn/(zn-zf)  1

//======================================================================///

inline void Matrix4x4OrthoOffCenterRH(CLMatrix4x4& mOut, float l, float r, float b, float t, float zn, float zf)
{
	mOut.m01 = mOut.m02 = mOut.m03 = mOut.m10 = mOut.m12 = mOut.m13 = mOut.m20 = mOut.m21 = mOut.m23 = 0.0f;
	mOut.m00 = 2.0f / (r - l);
	mOut.m11 = 2.0f / (t - b);
	mOut.m22 = 1.0f / (zn - zf);
	mOut.m30 = (l + r) / (l - r);
	mOut.m31 = (t + b) / (b - t);
	mOut.m32 = zn / (zn - zf);
	mOut.m33 = 1.0f;
}

//======================================================================///

//	      	2/(r-l)      0            0           0
//	      	0            2/(t-b)      0           0
//	      	0            0            1/(zf-zn)   0
//	      	(l+r)/(l-r)  (t+b)/(b-t)  zn/(zn-zf)  1

//======================================================================///

inline void Matrix4x4OrthoOffCenterLH(CLMatrix4x4& mOut, float l, float r, float b, float t, float zn, float zf)
{
	mOut.m01 = mOut.m02 = mOut.m03 = mOut.m10 = mOut.m12 = mOut.m13 = mOut.m20 = mOut.m21 = mOut.m23 = 0.0f;
	mOut.m00 = 2.0f / (r - l);
	mOut.m11 = 2.0f / (t - b);
	mOut.m22 = 1.0f / (zf - zn);
	mOut.m30 = (l + r) / (l - r);
	mOut.m31 = (t + b) / (b - t);
	mOut.m32 = zn / (zn - zf);
	mOut.m33 = 1.0f;
}

inline float GetZnFromProjectionMatrix(const CLMatrix4x4& mat)
{
	return -mat.m32 * mat.m23 / mat.m22;
}

inline float GetZfFromProjectionMatrix(const CLMatrix4x4& mat)
{
	return abs((mat.m33 - mat.m32) / (mat.m22 - mat.m23));
}

NAMESPACE_LABOR3D_END

#endif // __LABORMATRIX4X4_H__

